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In an effort to develop a chemically reactive interaction potential suitable for application to the study of con- 
ventional, organic explosives, we have modified the diatomic AB potential of Brenner et al. such that it exhibits 
improved detonation characteristics. In particular, equilibrium molecular dynamics (MD) calculations of the 
modified potential demonstrate that the detonation products have an essentially diatomic, rather than polymeric, 
composition and that the detonation Hugoniot has the classic, concave-upward form. Non-equilibrium MD cal- 
culations reveal the separation of scales between chemical and hydrodynamic effects essential to the Zel'dovitch, 
von Neumann, and Doring theory. 



I. INTRODUCTION 

For over 15 years, the Reactive Empirical Bond-Order 
(REBO) potential of Brenner et al. LLC, or variations thereof, 
have been used for Molecular Dynamics (MD) simulations of 
detonation. It has been shown to follow the Chapman- Jouguet 
(CJ) theory and even that of Zel'dovich, von Neumann, and 
Doring (ZND) [2]. Being an MD model, it has at least two 
major shortcomings, namely its spatial and temporal scales. 
As computers become more powerful and codes and algo- 
rithms become more advanced, some of the restrictions that 
limit these scales can be loosened, allowing for more realistic 
behavior to be modeled. 

While providing an atomic scale model of a reactive mate- 
rial, the empirical and classical model of Brenner et al. falls 
short of including all aspects of an accurate representation. 
Additionally, the spatial and temporal scales of real explosions 
make MD simulations a large computational task. But MD is 
still a useful tool for probing the characteristics of the detona- 
tion phenomenon, and the REBO potential is one of the best 
at balancing realism in the potential with accessibility to the 
large-scale; yet there is room for improvement. 

MD has certain conceptual advantages over hydrodynam- 
ics approaches that are parameterized to match the behavior of 
real high explosives. Even though the latter models can mimic 
real experiments on the proper spatial and temporal scales, 
they make assumptions about the reaction rate and multiphase 
equation of state (EOS) to do so. Some of these assumptions 
are parameter fittings, which may not elucidate any new phys- 
ical intuition about detonation of the high explosive (HE) in 
question. In comparison, MD simulations depend on the pa- 
rameterization of an interaction potential, which is arguably 
easier to connect directly to physical considerations than is a 
multiphase and multi-species EOS. 

REBO has been used by many groups to model a variety of 
parameterizations and experimental configurations. However, 
a major criticism of this potential is that it has a thin reaction 
zone (~ 100 A) relative to typical real high explosives (~ 
1 mm) fToh . In previous work HHH] several unconventional 
characteristics of the default parameterization of Brenner et 
al. [Q]] (called Modell as in 0]) are made evident: 

• Modell displays nearly instantaneous dissociation upon 
compression by an unsupported detonation. Its entire 



reaction zone is characterized by a dissociative state. 

• Its reaction zone is thin. 

• It readily allows for clustering. 

• Its CJ state is highly compressed. 

• It has clearly non-hyperbolic equilibrium Hugoniots 
(H) in P-v space. 

None of these characteristics are proven to be unrealistic. 
In fact for some primary high explosives, the plasma-like state 
seen in Modell at CJ may be realistic |5]. Clustering, partic- 
ularly of carbon, is a real phenomenon in the thermal decom- 
position of some HEs laUD- However, successful models of 
conventional explosives have either assumed or suggested that 
more molecular states with hyperbolic Hs in P-v space 
OHOjlEII] are typical and that compressions at CJ conditions 
should be « 25% (for example PETN fU) rather than « 43% 
as with Modell yfl. Given that Modell can model a dissoci- 
ated state, it would be useful to show that a more molecular 
state can also be modeled. The task of this work is to make 
modifications to the Modell potential in order to accomplish 
this goal. Although, for the last feature in the above list, it is 
difficult to predict what changes would adjust it, an attempt 
to address the former features is made by making adjustments 
based on physical reasoning (see Sec. HH). In Sec.[lVlthe ef- 
fects that the changes to the model of Brenner et al. have on 
the thermodynamic properties are investigated by the meth- 
ods outlined in Sec. [TTTJ Notable results of non-equilibrium 
MD simulations are studied in the companion paper. 



II. CHANGING REBO: PHYSICAL AND AESTHETIC 
REASONING 

In this section we motivate changes to the Modell poten- 
tial to form a new potential (called ModellV after the naming 
scheme of White et al. |0]), the total bonding energy of which 
takes the form 
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TABLE I: The components and parameters used in Eq. \\\ 
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The parameters for EqQ]can be found in TableU Before men- 
tioning the physical basis for the modifications made to Mod- 
ell, let us list the features that Modell already takes into ac- 
count. In tests using Modell, the parameterization is set to 
a valence of one, which should prefer dimers because one 
atom's sharing of electrons with one other fills both atoms' 
valence shells. The bonding energy takes the form 
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See Table I in the errata of Brenner et al. for the functions 
and components of Eq.O In Modell this bond-order behavior 
is implemented by reducing the strength of the attractive term 
of a Morse potential used to model such covalent bonding. 
This is accomplished with a bond-order coefficient {Bij ) that 
varies from one to zero with increasing proximity and number 
of neighbors. 

Rice et al. found that Modell allowed for trimer formation 
and subsequently strengthened this bond-order effect by in- 
creasing the repulsive term as well in order to better model 
bond saturation so that a particle became less likely to bond to 
more than one other |[l3l[lJ]. The contribution of Rice et al. is 
a simplified adjustment for what can be a complicated inter- 



action. For instance, the Bij coefficient does not take spin or 
degeneracy into account. However, this type of correction is 
similar to the "bond saturation" term that is used in the more 
sophisticated and calibrated ReaxFF potential of van Duin et 
al. 

In combination with this bond-order functionality, the vari- 
ation of bond well depths allows for reactive chemistry. If i 
and k are different types of atoms, their respective well depths 
with j will be different. If j's is deeper with k than with i, k 
will require less kinetic energy to displace i than i would have, 
had their rolls been reversed. That difference in energy will 
be converted into the kinetic energy of the system, an exother- 
mic reaction. This is the type of reaction modeled by Modell, 
2AB A 2 + B 2 + 2Q, where Q = D^ A/BB - D AB is the 
exothermicity of the reaction. The amount of energy needed 
to dissociate AB is D AB . To dissociate AA or BB requires 
£)AA/bb These num bers ignore the small contribution of the 
V v dw term m Eq.[2 

In Modell the van der Waals (vdW) interaction is repre- 
sented by a Lennard- Jones (LJ) form. It is connected with a 
terminating third-order polynomial with negative curvature in 
the short range domain so that the overly repulsive twelfth- 
order term does not compete with the Morse potential, which 
handles the covalent bonding and Pauli Exclusion repulsion. 
The LJ is parameterized to allow for solid lattice formation at 
low temperature. It turns out that the third-order inner spline 
introduces an artifact in Modell since it has a section of posi- 
tive slope, which allows it to trap particles (see Fig.Q]) because 
of the resulting attractive force. It is only first-derivative con- 
tinuous at the spline point and is cropped at the cutoff point. 

Modell simulates well the repulsion between individual 
atoms separated by less than their equilibrium distance — 
which depends on the number and proximity of neighbors to 
the pair. It does not, however, represent the electrostatic repul- 
sion between dimers that occurs as their charge clouds over- 
lap on approaching each other while still being too far away to 
rearrange bonds, ~ 3 A. As a result there is no significant in- 
teraction between molecules at the range and magnitude that 
gives rise to high pressure dense molecular HE product fluid 
mixtures 11611. T here are some good models for this, for ex- 
ample, ZBL fl7ll and ReaxFF 1151 . In an attempt to improve 
upon Modell and its results, we design a new model (Mode- 
ll V), which incorporates, along with the aforementioned con- 
tribution of Rice et al, the following modifications: 

The inner spline is replaced with a repulsive core, which 
smoothly connects the negative (the zero of the potential en- 
ergy is the dissociated state) LJ curve to a constant plateau in 
the region in which the Morse potential dominates. This inner 
core has no sections of positive slope — that is, sections of at- 
tractive force — and therefore cannot trap atoms, and it crudely 
models electrostatic repulsion. 

The reason for connecting to a plateau is that it allows one 
to dictate the depth of the bonding potential. It also simplifies 
the definition of a covalent bond. Any two particles that are 
within the defined bond distance (see the f c term in Table [J) 
with a radial component of kinetic energy less than the height 
of the plateau are considered bonded. This definition also dis- 
tinguishes bonded particles from particles that are merely con- 




FIG. 1: Total bond energy in a ID, three-particle, A-B-A, Modell 
interaction in which the distance between the B and one of the As is 
fixed at the bond's equilibrium distance vs. the distance between B 
and the remaining A (short dash). V vc lw for any two atoms (solid). 
Total bonding energy less the V v aw contributions (long dash). All 
curves are subtracted by their corresponding values at the position of 
the local minimum of the first curve for comparison purposes. This 
shows that the section of positive slope in the REBO V V dw inner 
spline can trap particles. Without it, there is no trap. 



fined by neighbors, and it can be determined at any instant. 
The height of the plateau, 1 eV, is chosen to be comparable to 
the thermal energy at the CJ state for Modell (30, an addition 
that should decrease the reaction cross section, desensitizing 
the HE to initiation and, perhaps, adjusting reaction time and, 
thus, the width of the reaction zone. Perhaps, a more physical 
function would continue to increase monotonically to a finite 
value at r = and have V v dw be particle-type dependent. 

All of the spline points in ModellV have, at least, contin- 
uous second derivatives, with the exception of the outermost 
cutoff point. To terminate the V v dw term m the long range, a 
Holian-Evans spline ifTill . which is second derivative continu- 
ous at its inner spline point but only first derivative continuous 
at its outer one, is used. A computational advantage for us- 
ing the Holian-Evans spline is that the reach of the ModellV 
potential is shorter (« 5.19 A) than Modell's (« 7.32 A). 
Fig. [2] shows a comparison of the two vdW potentials and 
Fig. [3] shows the corresponding forces. One can see that the 
ModellV V v dw term is still smooth in the force. As with 
the spline points in the V v dw term, the cutoff function in the 
bond-order function and Morse potential is also replaced so 
that it is second-order continuous. Second derivative conti- 
nuity helps energy conservation in the MD simulations (and 
allows longer time steps to be taken). 

Another difference between ModellV and Modell is the 
depth of the metastable covalent wells. With the addition of 
the repulsive core, ModellV fails to detonate with a well as 
deep as the default value for Modell. The metastable well is 
raised by 1 eV, such that D AB = 1.0 eV, D AA/BB = 5.0 eV, 
and therefore, Q — 4.0 eV. Comparing the potential energy 
surface for the linear, symmetric, metastable configuration for 
Modell (see Fig. [4]) to the one for ModellV (see Fig. [5]), we 
notice that the depression near (1.2 r e , 1.2 r e ) diminishes 
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FIG. 2: (Color online) van der Waals (V v dw) potential vs. the in- 
teratomic distance (r^). Spline points are indicated by changes in 
color-dashing. Modell V v dw is made of a third-order polynomial 
(cyan-solid) connected to a Lennard- Jones (LJ) form (mauve-short 
dash), which is cropped at a finite distance at which the slope and 
value of Vvdw are nearly zero. ModellV V V dw starts as a con- 
stant (red-long dash) and is connected with a fifth-order polyno- 
mial (green-long short dashes) so that the spline points are second- 
derivative smooth, to a LJ form (blue-short short long dashes) at 
zero, below which it coincides with the Modell LJ until the inflec- 
tion point, at which it is connected to a curve (black-long long short 
dashes) that is brought smoothly to zero within a finite distance, a 
Holian-Evans spline fl8ll . 
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FIG. 3: (Color online) Forces that correspond to the potentials in 
Fig.0 

somewhat for ModellV and that the barrier to be overcome 
for one particle to displace another via an end-on attack is sig- 
nificantly increased. All of these changes have notable effects 
on the EOS and detonation properties of the model HE. The 
final difference is that the masses of the particles are changed 
so that they are different from each other. The changes to 
Modell to form ModellV may be summarized as follows: 

• A bond-order coefficient is applied to the repulsive term 
in the Morse part of the interaction potential. 

• The inner spline of the V v dw term is replaced with a 
repulsive core, which comprises a plateau connected to 
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FIG. 4: Potential energy surface for the linear, symmetric, metastable 
configuration of Modell atoms. The contours are spaced every 
0.1 eV. ta b is the distance between and A and B atom, tb-a is 
the the distance between the same A atom and a different B atom. r e 
is the equilibrium distance for the metastable covalent interaction. 
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FIG. 5: Potential energy surface for the linear, symmetric, metastable 
configuration of Modell V atoms. The contours are spaced every 
0.1 eV. 



an LJ form with a monotonically decreasing fifth-order 
polynomial. 

• The order of all of the splines is increased such that all 
but one spline point (the outermost cutoff) have smooth 
second derivatives. 

• The masses of the atom types are changed so that they 
are no longer equal. 

• The depth of the binding energy for the metastable 
bonds is made shallower. 



III. METHODS 

Investigating the two REBO potentials of which this paper 
is a study, Modell and ModellV, we use SPaSM 3.0 fl9h to 
carry out the MD simulations, of which there are two main 
types, microcanonical ensembles (NVE) and non-equilibrium 
MD (NEMD). The NVE simulations are used for testing the 
thermochemical properties of a model and NEMD are used to 
investigate its detonation behavior. Both utilize the leapfrog 
Verlet method to advance the atoms' positions and momenta. 

We conduct two sub-types of NVE simulations. The first 
is meant to find the equilibrium Hugoniots and are similar 
to those described by others O, LLJ, 12011. The new model 
requires a shorter time step St « 0.16 fs than does that of 
Brenner et ah 01, St « 0.25 fs because of the relatively 
high curvature of the fifth-order spline in ModellV's V v dw 
term. Its zero-pressure configuration is also slightly differ- 
ent. The length of the sides of the rectangular unit lattice 
cell are l x = 6.19122 A and l z = 4.20538 A. The cell con- 
tains two dimers in a herringbone configuration. The angles 
that the dimers make with the horizontal are ±27.7109°. The 
atoms are placed 0.59976 A from the centers of their respec- 
tive dimers, which are positioned at (1/4 l z , 1/4 l x ) and (3/4 l z , 
3/4 l x ) from the lower left corner of the cell. These values are 
determined by isothermal-isobaric Monte Carlo simulations. 

The second type of NVE simulation, the cookoff, is used 
to determine the reaction rate. The initial conditions are the 
same as with the first type of NVE except that the constituents 
are lOOx 100 cells 2 of AB not 25x25 of AA and BB. The 
time step is also drastically reduced so that good statistics can 
be found for measurements taken during rapid reactions. The 
time step varies among cookoff simulations, 6 x 10 -4 fs < 
St < 0.02 fs. 

One of the purposes of creating a new potential is to expand 
the reaction zone by lowering the reaction cross section. This 
requires larger NEMD simulations to be performed. As Rice 
et al. have done fl3ll for computational efficiency, we try to re- 
duce the amount of pre- shocked material modeled by tacking 
initial state material onto the end of the sample as the shock 
front approaches. To reduce surface effects, we introduce at 
the surface a one-cell-thick layer of two new particles, C and 
D, that have all of the properties of A and B, respectively, ex- 
cept that they start frozen and their velocities are not updated. 
When the shock front nears the frozen layer, the frozen layer 
is converted to A and B and is thermalized along with the new 
material. Another frozen layer is tacked onto the end. The 
size of an NEMD simulation is discretely dependent on time. 
The largest simulation contained over three million atoms. 

Another effect of decreasing the reaction cross section is 
that the incubation time, before which detonation occurs, ex- 
tends. The sensitivity to impact is thus lessened. To over- 
come this, we initially overdrive the simulation with a piston 
that moves into the material at a velocity u pstn « 4.9 km/s. 
After a detonation front has seemingly been established, it is 
smoothly backed off by following a sine- shaped deceleration 
over 200 time steps to a desired velocity. The time step is the 
same as with the first type of NVE simulation, St « 0.16 fs. 
Most of the results of the NEMD simulations are reported in 
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FIG. 6: Snapshot of a section of a shock front for an under- supported 
detonation using the ModellV potential. Shock is propagating up- 
ward. Particles are shaded by atom type. 



the companion paper. 



IV. COMPARISON OF THERMODYNAMIC PROPERTIES 

We start the comparison of the properties of Modell and 
ModellV by analyzing how the ModellV material behaves 
when compressed by the passage of an under-driven detona- 
tion front. From the snapshot of an NEMD simulation using 
ModellV (see Fig. 0), one can see on the left half that many 
unreacted dimers are lining up horizontally. This is represen- 
tative of uniaxial compression. Farther back it has melted, and 
farther back still reacted products are evident. On the right 
half, the dimers do not seem to line up before reacting. It is 
made evident in the companion article that this is an effect of 
detonation instability and the propagation of transverse waves. 
When comparing this to a snapshot of detonation in Modell 
(see Fig. 1 in the paper of Heim et al. OD), one notices that, 
even at the right side of the sample, ModellV seems to hold its 
molecular identity better than Modell during compression by 
a detonation front that is not overdriven. It seemingly displays 
a greater resilience to dissociation. Both the right half of the 
ModellV snapshot and the Modell snapshot lack a significant 
induction zone. 

To better compare the widths of the detonation fronts of the 
two models, profiles of the average z -component of particle 
velocity from a critically supported detonation of ModellV 
are plotted. Figure [7] indicates that the reaction zone is about 
700 A wide since that is the position at which the curves settle 
down to the constant u P j . With the same method, the width 
of the reaction zone for Modell was determined to be about 
300 AH. 

Although ModellV's reaction zone is still small compared 
to that of real explosives, it can be widened. To widen the 
reaction zone farther, one can create a reaction that requires 
more steps. One also can have the reaction increase the num- 
ber of moles of material. The resulting expansion in volume 
would be another driver for the shock wave yfl. Yet another 
effect of this would be to introduce an entropic penalty to 
back-reaction and make the model more closely approximate 
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FIG. 7: Overlap of profiles of the z-component of the particle ve- 
locity for a critically supported NEMD simulation using ModellV 
parameterized such that = 1.0 eV and Q — 4.0 eV. A constant 
line is drawn at the CJ value of the particle velocity (u P j), at which 
the driving piston is moving. Using a critically supported NEMD as 
opposed to an unsupported one drastically reduces the transient time 
from initiation to steady state [3fl. 



TABLE II: Determined thermodynamic values at CJ for ModellV, 
where v is the specific volume, u s is the shock velocity, u p is the 
particle velocity, T is the temperature, U is the potential energy, E is 
the internal energy, P is the pressure, and A is the degree of reaction. 
Subscript j indicates the CJ value, Subscript indicates the value at 
the initial state. Averages are per particle. 



Vj/VQ 


0.789851(39) 


U s j 


13.87868(50) 


U p j 


2.91658(44) 


(k B Tj) (eV) 


0.9185(84) 


(Uj) (eV) 


-0.8631(72) 


(Ej) (eV) 


0.0554(12) 


Pj (eV/A 2 ) 


0.83846(16) 


A; 


0.8548516(36) 



the ZND assumption of irreversibility. 

To find the critical piston velocity in the preceding analysis, 
the CJ state for ModellV is found by conducting NVE simula- 
tions and seeking the values of E and v that satisfy the Hugo- 
niot jump conditions O, [lj, [20|] . For Modell Vj/vo « 0.57 
and u S j « 9.7 km/s |3]. From the results (see Table [II]), one 
should notice that for ModellV the CJ state is less compressed 
than for Modell and that u S j is faster. Compared to con- 
ventional explosives, which typically have specific volumes 
Vj/vo « 0.75 and shock velocities u S j « 6 km/s (for ex- 
ample PETN tun, ModellV detonation fronts are over twice 
as fast and its CJ state is slightly less compressed. Note, how- 
ever, that it has been shown for these REBO potentials that 3D 
simulations have lower, and therefor more realistic, velocities 
than do 2D ones 0,|2l|] . Heim et al. JH] show how closely the 
conditions of the final state in a detonation of a Modell mate- 
rial match those of the CJ conditions. For contrast in Fig. [5] 




FIG. 8: Hugoniots for Modell (circles) and ModellV (squares). Solid 
lines are a guide to the eye for Modell and a fit for ModellV. Rayleigh 
lines are determined by the initial conditions and slopes. Dotted lines 
represent Rayleigh lines the slope of which are determined by the CJ 
value of the detonation velocity. Slopes of dash-dotted lines are de- 
termined by the average velocity of the shock waves in unsupported 
detonations. Boxes are magnifications. 



one can see that the Rayleigh line for the unsupported simula- 
tion is significantly steeper than for that of the CJ state. 

In ZND theory this would require that the final state for 
ModellV is either at a strong point or a weak point. The strong 
point is unstable and will only be a solution to the conserva- 
tion equations if the detonation is overdriven. In the strict 
ZND theory there is no path to the weak point. ZND is, how- 
ever, based on assumptions not realized by ModellV. Mode- 
llV has a reversible reaction, the pathway of which includes an 
endothermic dissociation. Endothermic steps can be responsi- 
ble for weak point final states |2]. Perhaps, the repulsive core 
introduced in ModellV, by sheltering the dimers, prevents the 
reduction of the activation energy of the dissociative step and, 
hence, increases that step's reaction time, thus making it a 
stronger contributer to the overall reaction rate and causing 
the final state to be at a weak point. As is shown in the ac- 
companying paper and can be inferred from Fig. [6J however, 
a ID theory is not totally appropriate for this model. 

In order to confirm that the modifications to the REBO 
model in the ModellV potential do, indeed, maintain a dimer- 
ized state, the radial distribution function (RDF) at the CJ state 
is plotted. In Fig. [9] the CJ state of ModellV is compared to 
that of Modell. From the RDF for ModellV, one can see that, 
after the peak at r = r e , the curve drops to nearly zero while 
the analogous region in the RDF of Modell is closer to unity. 
This is indicative of a molecular state for ModellV and a dis- 
sociative state for Modell. From Fig. [TUJ it is clear that the 
snapshot of the ModellV CJ state supports the finding of the 
corresponding RDF. The snapshot for Modell shows many 
more clusters and dissociated atoms. 

For these images, a bond is defined for each potential. Two 
particles, i and j, are considered bonded if E^^j + KEu^j < 
E c and < r c , where E^j is as in Eq.[T]for ModellV and 
Eq.[2]for Modell except that the outer sum is over only i and j. 



0-0 A-A, ModellV 
o-o A-B, ModellV 
Q- £] A-A, Modell 
a-a A-B, Modell 




FIG. 9: Radial distribution functions for the CJ states of ModellV 
and of Modell. 
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FIG. 10: Snapshots of the CJ state for ModellV (upper) and Modell 
(lower). Particles are marked by bond type. Gray atoms are in un- 
reacted dimers; striped are unbonded, dissociated atoms; black are 
clustered atoms; and white are reacted dimers. 



KE^ij = - fi \ (vi — Vj) • rij/rij\ . [i is the reduced mass of 
i and j. For ModellV E c = ec, and r c = 72 (see Table IJ). For 
Modell E c is the peak of the inner spline on the V v dw term 
and r c is the location of that peak (see Fig.O. It can be shown 
that for Modell there exists no minimum above and within 
these cutoff values for all values of B^. To be considered 
part of a dimer, a particle must be bonded to only one other 
particle that is bonded to no other. If the particle types are the 
same, it is a product dimer; if different, a reactant dimer. To be 
part of a cluster, a particle must be bonded to either multiple 
particles or to one other that is bonded to multiple particles. 
To be dissociated, a particle must be bonded to no other. 

With a bond defined, cookoffs that will reveal the reaction 
rate and the activation energy (E a ) for each model can be sim- 
ulated. This is done for two reasons, to compare ModellV to 
continuum models and to show another fortuitous difference 
between Modell and ModellV. The cookoffs are NVE simu- 
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FIG. 1 1 : Reaction rate of ModellV for various internal energies and 
volumes vs. temperature. Lines are fits to the Arrhenius form for con- 
stant volume. Their slopes are the negative of the activation energy. 
t u is a unit of time and equals 10.180505 fs. 



FIG. 12: Reaction rate of Modell for various internal energies and 
volumes vs. temperature. Lines are fits to the Arrhenius form for 
constant volume. Slope is the negative of the activation energy. t u is 
a unit of time. 



lations started at different values of compression and internal 
energy but all with the initial AB chemical composition. As a 
simulation is started, thermal energy is partitioned among the 
different modes. After this has occurred and the reaction has 
progressed somewhat, the reaction rate is sampled over a short 
time and the temporal average of the rate and its standard er- 
ror are recorded for the current simulation. In Fig. [TT] the data 
are plotted for a series of simulations using the ModellV po- 
tential. For each value of v, the data are fit to an Arrhenius 
reaction rate of the form Il25n 



(3) 



A = (1 — A)Aexp 



where A is the degree of reaction defined computationally as 
the ratio of reacted (defined as above) material to total amount 
of material. A is the reaction rate and A is the frequency factor. 
The analogous plot for Modell can be found in Fig.fT2l 

When contrasting the two figures, one notices that, over a 
similar range of volumes, ModellV maintains a better fit to 
the Arrhenius form than Modell and that the calculation of 
E a remains relatively constant when compared to the E a of 
Modell. A probable contribution to the former observation is 
the manner in which a reaction is defined. At compressions 
typical in detonations, Modell has a larger number of atoms 
in clusters and free states, which are not counted as reacted. 
Only stable dimers are. One may consider a cluster of two 
A atoms tightly bonded together with a third, B-type particle 
loosely bonded to one of them as a state more reacted than one 
in which A is tightly bonded to B and loosely bonded to the 
other A. However, by our definition of reaction, all particles 
in a cluster are considered unreacted and only contribute to 
the denominator of the calculation of A. 

E a for Modell decreases with decreasing v. This is be- 
cause at higher compressions an atom can have more neigh- 
bors within the bonding distance than with ModellV The 
bond-order coefficient then lowers the attraction between it 



and its neighbors, making it easier to break apart any exist- 
ing bonds. For ModellV, the repulsive core keeps neighbors 
away from dimers so that the bond-order coefficient and, thus, 
the dissociation energy remain unaffected by the same level 
of compression that would otherwise cause Modell's dissoci- 
ation energy to be reduced. For ModellV there is a higher ten- 
dency than for Modell for the compression to use up the space 
between dimers than that between the dimers' constituents. 
This is because for the former case a third particle would have 
to climb the repulsive core that was added to the V v dw term 
before it got within the range of the cutoff function f c , where 
it can diminish the bond-order coefficient . In the latter 
case a third particle would only need to overcome the barrier 
of Modell's inner spline (see Fig. [2]), which is small compared 
to the temperature for most situations considered in reaction 
and detonation. For ModellV the activation energy goes up 
slightly with decreasing volume. The exact reason for this is 
not known at present, but simulations with ReaxFF on the high 
explosive RDX also show this trend [22]. E a for ModellV is 
roughly 1 eV greater than for Modell. This is commensurate 
with the height of the repulsive core. 

Even though we show that ModellV fits an Arrhenius reac- 
tion rate in NVE cookoff simulations, we would like to know 
that this is the type of reaction that takes place in the reac- 
tion zone. We can use that reaction rate along with the von 
Neumann spike (vNs) state to predict a A profile. Since we 
do not have the full EOS throughout the reaction zone, we 
make some assumptions that will tend to shorten the width of 
the profile. Let's assume that u p increases linearly and that 
v and P are constant such that T(A). We conduct a series 
of shock tube simulations of a substance similar to ModellV 
except that Q is set to zero. We find the value of piston ve- 
locity that produces the constant zone the state of which falls 
on the theoretical 1Z from Fig. [8j This is the vNs state. The 
value of the specific volume at the vNs state is v n w 0.558 vo 
and temperature T n « 0.25//^. From FigQj] we estimate 
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FIG. 13: Theoretical profiles for the degree of reaction (A), temper- 
ature (T), and reaction rate (A) for ModellV, assuming a constant 
pressure and density, max [A] =8.94 THz. 

the reaction rate parameters (see Eq. [3]) at the vNs state to be 
A n « 3.6 and E an ^2.0. Inserting the following equation 
for temperature 

T(A)=T n + (T j -T n )A (4) 

A 3 

into Eq.[3]and using the change of variable 

dX dX dt • / u V j\~ 1 

Tz = mT z =x \ u --—) ' (5) 

we generate a profile by iterating backward in space via a third 
order Runge-Kutta method (see Fig [13]). The variables with 
the subscript j are the CJ values listed in Table HH 

The theoretical A profile reaches Xj at z « —70 A, much 
closer to the front than the —700 A that we measure for the 
CJ value of the particle velocity in Fig. [71 This does not prove 
that the reaction rate is Arrhenius in the reaction zone, but it 
does show that the order of magnitude of the measured reac- 
tion zone width is significant compared to the lower bound of 
the theoretical width. Had it been otherwise, this would have 
shown the reaction rate to not be Arrhenius and thus a differ- 
ent form, for example, a flame front. 

Replacing the assumptions in the analysis with more ac- 
curate ones, for example, variable v, would tend to lengthen 
the theoretical width because the expansion of the material 
would tend to cool it and slow the reaction rate. In lieu of the 
full EOS, we could have used data from an NEMD to fill in 
the state within the reaction zone. However, we would only 
be testing internal consistency. This comparison also suffers 
from an assumption of steady-state one-dimensionality. We 
have seen that the NEMD simulations are 2D. As we show in 
the companion paper, transverse waves traverse the reaction 
zone. This causes the profile to be time and x dependent. The 



EOS throughout the reaction zone is currently being calcu- 
lated for normal mode stability analysis II 2411. which assumes 
an Arrhenius reaction rate. 

V. CONCLUSION 

In this paper it was shown that the changes to the Mod- 
ell version of the REBO potential have reduced the amount of 
clustered and dissociated atoms at the CJ state (by RDF) if not 
throughout the reaction zone (by a snapshot of an NEMD sim- 
ulation). The CJ state of ModellV is at a less compressed state 
and ModellV's reaction zone is wider than Modell's is. Mod- 
ell's unsupported detonation velocity is much better predicted 
by CJ theory, where ModellV's propagates significantly faster. 
In the companion paper we investigate ModellV's relative dis- 
parity more closely. The cookoff simulations reveal that Mod- 
ellV fits an Arrhenius reaction rate through a wider domain of 
compressions than does Modell and it has a higher and less 
density-dependent activation energy. 

It was put forth that this higher E a causes the final state 
to be at a weak point since it might make more significant in 
the reaction rate the endothermic step of dissociation. Given, 
however, the 2D shape of the front in Fig.O such a ID expla- 
nation is probably insufficient. Higher E a s should increase 
the likelihood of 1 and 2D instabilities in detonation waves 
HH [13]. The current value of w 2 eV is consistent with the 
presumed activation energies for many conventional HEs. 

The goal of the current research is the improvement of 
REBO, which we believe we have achieved, mainly, by in- 
troducing increased atomic repulsion, which reduces the reac- 
tion cross section. It was shown that ModellV's reaction zone 
is wider than with Modell and that the CJ state is a molecu- 
lar one. REBO's insensitivity to initiation is increased, and a 
thicker induction zone is introduced. The reaction rate now 
behaves in a more Arrhenius manner, and the product EOS 
behaves more like a polytropic gas as indicated by the hyper- 
bolic shape of its equilibrium H. Some measurements, for ex- 
ample, the adiabatic 7 — as is shown in the companion paper — 
and the amount of compression at CJ, are more commensurate 
with conventional explosives [10] . However, some are not too 
close, for example, the reaction temperature and shock veloc- 
ity. 3D simulations are needed to investigate how the changes 
truly compare to real experiments. 
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